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ABSTRACT 

Thick rubber components are employed by the Army to carry large loads. In tanks, rubber 
covers road wheels and track systems to protect roadways. It is difficult for design engineers to 
simulate the details of the hysteretic heating for large strain viscoelastic deformations. In this 
study, an approximation to the viscoelastic energy dissipated per unit time is investigated for use 
in estimating mechanically induced viscoelastic heating. Coupled thermo -mechanical 
simulations of large cyclic deformations of rubber cylinders are presented. The cylinders are 
first compressed axially and then cyclically loaded about the compressed state. Details of the 
algorithm and some computational issues are discussed. The coupled analyses are conducted for 
tall and short rubber cylinders both with and without imbedded metal disks. 


INTRODUCTION 

Rubber is employed to carry large loads in tires, gaskets, and tank track pads. It is also used to 
provide damping and system stability in complex mechanical systems such as helicopter rotors. 
In these applications, the rubber is typically stiffened by the addition of carbon black. The filled 
rubber tends to be a poor conductor of heat, yet it also exhibits very large hysteretic energy loss 
during cyclic loading. Also, the mechanical properties of rubber are strongly dependent on 
temperature. Faced with the above issues, designers interested in modeling the detailed response 
of complex-shaped rubber components need to be able to compute the coupled thermo- 
mechanical behavior of rubber. 

An example of the importance of understanding the thermo -mechanical response of filled rubber 
is given in a series of papers presented at the "Thirty Second Sagamore Army Materials Research 
Conference" held a Lake Luzerne, NY in 1985.'’^ In these papers, hysteretic heating, thermo- 
mechanical degradation, and fatigue of rubber-coated road wheels and tank track pads are all 
discussed. Uncoupled thermo -mechanical finite element analyses, and sensitivity studies were 
conducted with finite element codes. It was observed that the viscoelastic properties and the 
shape of the rubber solid are the most important factors in determining temperature rise.' The 
degradation studies indicate that the failure of cyclically loaded "rubber-like" polyurethane 
blocks depends on the hard segment transition temperature.^ Experiments were conducted which 
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proved that the large strain hysteretic heating rate does not correlate with the heating rates 
predicted using the popular complex modulus material data.^ Also determined is the fact that 
failure under cyclic loading can be "significantly different from that obtained in constant strain 
rate testing.""^ These conclusions suggest that detailed computational simulations of large strain 
dynamic loading of rubber-like solids, performed as part of a material degradation study, require 
accurate modeling of the large strain viscoelastic properties and a coupling of the mechanical and 
thermal models. 

In this paper, large deformation rubber viscoelasticity and heat transfer finite element models are 
coupled and employed to simulate the hysteretic heating of dynamically loaded rubber cylinders. 
The purpose of the work is to establish and test a computational tool for analyzing hysteretic 
heating in rubber components. 

The formulations for large strain rubber viscoelasticity and heat transfer employed in ABAQUS 
are outlined below to facilitate the description of the thermo -mechanical coupling performed in 
this study. Detailed information on these two formulations and their finite element 
implementation is available in the ABAQUS Theory Manual.^ Additional information on 
formulations for rubber viscoelasticity is available in the literature. Only moderate 
temperature changes were simulated in this study, so time-temperature superposition is not 
discussed. 


RUBBER ELASTICITY 

The coordinates of a material point in the reference configuration are indicated by 
x = {x, and in the deformed configuration by x = {x,,X 2 >-t; 3 }. The vector function 

X = x(X, t) defines the mapping of points from the reference configuration into the deformed 
configuration. The Lagrangian method is employed to describe the deformed solid's energy. 
The matematical tools employed are listed in the APPENDIX. The virtual work statement, 
excluding inertial effects, for a deformed solid of volume V and surface area S is 

SWj=js :5D^dV = \dv^tdS + \5v^ fdV (1) 

V S V 

where 5vis a virtual displacement vector, dWf is the internal energy due to the virtual 
displacement, is the virtual rate of deformation, t is the traction stress prescribed over S , 
and f is the body force vector acting within the volume V , and Vq is the volume in the 

undeformed state. When the solid is rubber, the internal energy in Equation (1) is expressed as 
follows. 



(dU - dU 
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where U is the strain energy density function defined in this effort as 
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where Gq is the instantaneous shear modulus and kg is the instantaneous bulk modulus. 

RUBBER VISCOELASTICITY 

The present analysis employs the large strain viscoelasticity model adopted in ABAQUS.^ The 
model assumes that stress relaxation can be described by a Prony series. The shear and bulk 
moduli are expressed in the form 

G{t) = G^+f^G,e^ and K{t) = +f^K, (4) 

i=l /=1 

where Gq = G{t = O), K^ = K{t = O), and denote relaxation time constants. The constitutive 

model employed is a BKZ-like model in which the deviatoric, s and hydrostatic, s ^ , stresses 
are determined as follows. 


s^(t) = s^(t)-5TM 


!=1 "f,- ^ 


and 




( 5 ) 

( 6 ) 


where g,. = G, / Gq , k, = A,. / , and F, (/ -^ ) is the deformation gradient of the deformed shape 

at time relative to the deformed shape at time t. Approximations are made when 

Equations (5) and (6) are integrated so that the entire history does not have to be convoluted. 
The approximations depend on the material constants being of Prony series form. Throughout 
the remainder of this paper only one term in the Prony series will be employed and the Prony 
series subscripts will be dropped. 


HEAT TRANSFER 

The variational statement of the energy balance equation for heat transfer, together with Fourier's 
law, for a deformed solid of volume V , and surface area S is: 

\p^^56dV + \^^yi — dV = \5erdV + \5eqdS (7) 

J/ dt dx dx y ^ 

where 6 is the temperature, Ug is the internal thermal energy, p is the mass density, q is the 
heat flux per unit area, r is the heat supplied per unit volume, k is a conductivity matrix, and 
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56 is a virtual temperature field satisfying the essential boundary conditions. Equation (7) is 
used to build the transient heat transfer finite element equations for the rubber solid. 

Finite element discretization of the mechanical and thermal variational statements, Equations (1) 
and (7), results in systems of time dependent ditferential equations which approximate the 
mechanical and thermal response of the solid. It is difficult to compute the viscoelastic energy 
dissipated per unit time for the constitutive model described above. This mechanically generated 
energy dissipation is required if one needs to couple the mechanical and thermal analyses. 
Below we investigate approximating the energy dissipation with a simple formula involving the 
viscous components of stress. The formula is of the same form as the dissipation formula for a 
Maxwell model. 


ANALYTICAL 

The mathematical relation for computing the energy dissipation per unit time for a constitutive 
model in history integral form involves a double convolution integral in time." To avoid such a 
difficult computation we exploit the fact that when the relaxation modulus is in Prony series 
form, the convolution integral in the BKZ-like constitutive model can be approximated with a 
finite difference equation that is similar in form to the finite difference equation for a Maxwell 
solid. This suggests that we explore the use of the Maxwell solid's dissipation function (which is 
easily computed) for approximating mechanically induced heating in a BKZ-like viscoelastic 
solid. 


FINITE DIFFERENCE EQUATIONS FOR VISCOUS STRESS COMPONENTS 

A simplified one-dimensional version of the finite difference equations (across a time interval 
= for tho evolution of the viscous stress, cr'', that is valid for the BKZ-like model 

described above is given by^ 
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where (7^ is a viscous stress increment computed at time using the relative deformation 
gradient between configurations at times and , and using the instantaneous modulus, Gq , 
of the material. The analogous finite difference equation for the viscous component of stress in a 
Maxwell solid (obtained by employing the trapezoidal method) is given by 
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where (c7°^.j-(7°) is a stress increment computed using the total strain (i.e., current 
configuration's strain) and using the instantaneous modulus of the material. 

For a time interval in which the stress increment terms, , in Equations (8) and (9) are 

approximately equal, or in which the stress increment terms have small values relative to the 

-At 

stress decay term, e ^ ol, the Maxwell model and the BKZ-like model will have stress 
histories that are approximately equal. This implies that the dissipation functions will also be 
approximately equal. 


DISSIPATION FUNCTION 


A deformation suddenly imposed at time / = 0 and which is subsequently maintained constant, 
produces stress relaxation. When the rubber's relaxation shear modulus is of Prony series form 
the energy dissipation per unit volume, r{t), during relaxation is of the form 


du: -t/;Uoi 

dt T 


( 10 ) 


where C/j'(t = 0'") is the viscoelastic shear energy density which results from the sudden 
deformation. Note, the rubber solid is incompressible in this study and there is no dissipation 
from the volumetric viscous stresses. In the case of a Maxwell solid, the dissipation per unit 
volume, during a time dependent process is of the form 




-2u:{t) 

T 


( 11 ) 


where UJ (t) is the viscoelastic shear energy density in the Maxwell solid's spring-dashpot leg. 
Since the finite difference stress evolution equations in the finite element code used here 
(ABAQUS) are similar in form to Maxwell solid equations, we selected to employ the 
dissipation function in Equation (11), for estimating hysteretic heating from large strain 
viscoelastic deformations. Values of were approximated by assuming the current value of 
the viscous stress, s”(t), is derived from a linear internal solid in a Maxwell link which has a 
time constant given by T . The resulting dissipation function is given by 


r(t) 


2gG,T 


> 0 


where g Gq is the viscous shear modulus. 


( 12 ) 
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FINITE ELEMENT EMPLOYED 


The ABAQUS finite element code was employed to perform the calculations. The code allowed 
us to run the mechanical and thermal algorithms simultaneously in time and to exchange data 
between the algorithms. We present the details for the 2D axisymmetric model (the procedure is 
similar for the 3D solid model.) 

The eight-node hybrid axisymmetric CAX8RHT element was used, see Figure I. The 
interpolation is biquadratic for the displacement field, and bilinear for the temperature and 
pressure fields. Reduced integration is also employed. The nodal variables employed in the 
thermal element to describe the interpolation are the temperatures, 0 , at the four comer nodes. 
The stress element has radial, , and axial, ii^ , displacements at the eight nodes and hydrostatic 
pressure, , variables at each of the four integration points. The heating rates, r{t), are applied 
at the four integration points. 


HYSTERETIC HEATING IN RUBBER CYLINDERS 

To investigate the application of the procedure described above, four cylinders were analyzed for 
hysteretic heating. Only moderate temperature changes were simulated so the elastic material 
constants were not treated as temperature dependent. Having the elastic constants independent 
of temperature simplified the calculations presented in this effort, but it is not a limitation of the 
algorithm. 


CYLINDER DIMENSIONS AND MATERIAL PROPERTIES 

There were two groups of two cylinders each. One group consisted of uniform cylinders and the 
other group had steel disks at their centers, see Figure 2. All cylinders had the radius, 
R = 0.0282 m . There were two cylinder heights in each group. The heights were 
H = 0.05 m, and 0.0125 m respectively. The cylinders were compressed between steel fixtures. 
The model simulates the case when a lubricant maintains the fixture-rubber interface as 
fi-ictionless. The internal steel disks were completely attached (bonded) to the rubber. The 
height and radius of the disks were 0.0025 m , and 0.0I4I m respectively. 

The rubber energy density was modeled as a Neo-Hookean solid (Equation (3)). The viscous 
behavior was described with one Prony term (Equation (4)). The material constants employed 
are representative of a filled rubber and are listed below. The viscoelastic constants for the 
rubber are Gq =2.310 MPa, kg = k^ = 2.000 * 1 0^ MPa , g = 0.3, k = 0.0 MPa, and T = 0. 1 s . 
The elastic constants for the steel are Young's modulus £’ = 206.8 G Pa and Poisson's Ratio 
V =0.3. The film heat transfer coefficients for the rubber-air and rubber-steel interfaces are 

= 5.44284Y/(°C m s) and =20934 //(°C w s), respectively. The remaining thermal 
material constants are shown in Table I. 
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PRESCRIBED AXIAL DISPLACEMENT AND RESULTS 


The deformation of the tall { H = 0.05 m ) uniform cylinder subjected to an axial displacement is 
shown in Figure 3. Since the boundary conditions are symmetric, only the half height (r,z)- 
quadrant is shown. The mesh refinement near the top and outer edges is to accommodate the 
heat transfer gradients at those locations. Each cylinder was subjected to an axial enforced 
displacement that produced large strain hysteresis. The half height enforced axial displacement 
for the tall cylinder is described as follows. The top end of the cylinder was ramped to a 
displacement of - 0.0045 m in 0.05 s . The top end of the cylinder was then forced to follow the 
prescribed axial displacement given by 

uXt) = -0.0045 - 0.003 sin(40.84/) m (13) 

The half height enforced axial displacement for the short cylinder is described as follows. The 
top end of the cylinder was ramped to a displacement of - 0.001 125 m in 0.05 6' . The top end of 
the cylinder was then forced to follow the prescribed axial displacement given by 

uXt) = -0.001 125 - 0.00075 sin(40.84t) m (14) 

where t is the time in seconds. The loading and the displacement of the top of the cylinder are 
shown in Figure 4 for a time interval of 1.0 s. As expected, the displacement curve demonstrates 
viscoelastic softening. A number of analyses were performed to test the nonlinear elastic, the 
viscoelastic, and the transient heat transfer finite element algorithms separately. Hand 
calculations and finite difference calculations verified that the algorithms functioned correctly. 
The thermal boundary conditions did not significantly affect the temperature fields computed. 
The following analyses were performed, employing the dissipation function described above, to 
investigate the non-uniform heating of rubber cylinders undergoing large strain dynamic 
deformations. 

Uniform Cylinders. The cyclic deformations given by Equations (13) and (14) were applied to 
the tall and short uniform cylinders indicated by the dimensions given above. The temperatures, 
at the center of the cylinders, (r = 0,z = O), as a function of time for the first 20 s of loading are 
shown in Figure 5. The frictionless rubber- fixture boundary condition allows the strains to be 
uniform in the cylinder. Rubber is a poor conductor of heat and the fact that the cylinders heated 
nearly uniformly regardless of height was expected. 

Cylinders with Imbedded Disks. It is difficult to estimate hysteretic heating in rubber solids of 
complex shape because coupled thermo-mechanical analyses are needed in regions of high strain 
gradients. Tall and short cylinders with imbedded steel disks were cyclically loaded to simulate 
the distribution of viscoelastic heating in a rubber solid when the deformations produce high 
strain gradients. The results obtained appear reasonable. Figure 6 shows the meshes on the 
reference and deformed shapes for the tall cylinder. The temperature distribution in the tall 
block after 20 s of dynamic loading is shown in Figure 7. Temperatures at points A, B, C, and D 
in Figure 7 are plotted as a function of time in Figure 8. The outer radial end of the internal disk 
(point C) is predicted to heat much faster than the rest of the cylinder. 
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Similar results were obtained for the short cylinders. However, there were some convergence 
problems that were overcome by re-meshing near the top outer corner of the steel insert and 
employing lower order displacement interpolations (the four node displacement element was 
employed.) The temperature distribution in the short cylinder is shown in Figure 9. Again, 
Figure 10 shows that the region of high strain gradient, located near the outer radial end of the 
internal disk (point C), has the most rapid rise in temperature. 


SUMMARY 

Accurate predictions of the strain and temperature distributions in rubber components, employed 
in dynamically loaded structures, are required to perform degradation studies. A procedure that 
couples a viscoelastic large dformation stress analysis with a heat transfer analysis was described 
with the use of the ABAQUS finite element code. A user subroutine was written to approximate 
the time rate of energy dissipation per unit volume and to pass this data from the stress routine to 
the heat transfer routine. 

The thermo -mechanical heating of tall and short uniform rubber cylinders (without internal 
disks) was computed using the coupling procedure. The viscoelastic material properties 
employed are valid for large strain deformations of rubber. The time dependent strains in the 
cylinders were uniform, and uniform heating was computed. The thermo -mechanical heating of 
tall and short rubber cylinders, containing internal steel disks, was also computed. The internal 
steel disks provided high strain gradients within the rubber cylinders and nonuniform hysteretic 
heating was observed. The analyses performed suggest that the coupling procedure should be 
considered for further development as a design tool for rubber degradation studies. 
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APPENDIX 


The following notation is included to assist the reader with the description of the ABAQUS finite 
element algorithm for rubber viscoelasticity. 


A,. (i = 1,2,3) 
Xi(i= 1,2,3) 

X = x(X,t) 


coordinates of a material point in the reference configuration, 
coordinates of a material point in the deformed configuration, 
vector mapping between the reference and deformed configurations. 


F = ^ 

ax 


deformation gradient. 


J = det(F) 

F = / 

b = fF 


/><Kb) 

^2=ife)-»-(55)) 


determinate of F which measures volume change. 

deformation gradient scaled for volume change. 

left Cauchy Green strain tensor. 

first strain invariant (adjusted for volume). 

second strain invariant (adjusted for volume). 


5u 


5L = 


ddu 

dx 


displacement, 
gradient of displacement. 


5D = -(5L + 5L^) 
A:B 

& = 5D--5e''"^I 
3 

1, 

p = — 1 : s 
3 


rate of deformation. 

rate of deformation computed using virtual displacement 5\ . 
scalar product of matrices A and B . 
volumetric strain rate. 

deviatoric strain rate, 
pressure stress (hydrostatic). 
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Table I. 


Thermal properties. 


Property 


Rubber 

Steel 

Conductivity, J/(°Cms) 

K = 

0.20934 

45.83379 

Density, kg/m^ 

P = 

1000. 

7849. 

Specific heat, J l{kg °C ) 

c = 

2093.4 

460. 

Expansion, (°c) ' 

a = 

SO.^IO"® 

12.*10"® 
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FIGURE CAPTIONS 


Figure 1. Axisymmetric CAX8RHT element. 

Figure 2. Rubber cylinder, finite element mesh, fixture, and steel disk. 

Figure 3. Deformation of a tall uniform cylinder (FI = 0.05 m). 

Figure 4. Viscoelastic softening of a tall uniform cylinder (H = 0.05 m). 

Figure 5. Increase in temperature as a function of time at the center of each uniform cylinder. 

Figure 6. Deformation of a tall cylinder with an internal disk (H = 0.05 m). 

Figure 7. Temperature distribution in a tall cylinder with an internal disk (H = 0.05 m). 

Figure 8. Temperature as a function of time for a tall cylinder with an internal disk (H = 0.05m). 

Figure 9. Temperature distribution in a short cylinder with an internal disk (FI = 0.0125 m). 

Figure 10. Temperature as a function of time for a short cylinder with 
an internal disk (H = 0.0125 m). 
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Figure 1. Axisymmetric CAX8RHT element. 
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Figure 3. Deformation of a tall uniform cylinder (H = 0.05 m). 
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Figure 4. Viscoelastic softening of a tall uniform cylinder (H = 0.05 m). 
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Figure 5. Increase in temperature as a function of time at the center of each 

uniform cylinder. 



Figure 6. Deformation of a tall cylinder with an internal disk (H = 0.05 m). 
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Steel disk 

Figure 7. Temperature distribution in a tall cylinder with an internal 

disk (H = 0.05 m). 



Time (s) 

Figure 8. Temperature as a function of time for a taU cylinder with an internal 

disk (H = 0.05m). 




Figure 9. Temperature distribution in a short cylinder with an internal 

disk (H = 0.0125 m). 


Figure 10. Temperature as a function of time for a short cylinder 
with an internal disk (H = 0.0125m 





